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We identify and analyze quasiperiodic and chaotic motion patterns in tlie time evolution of a clas- 
sical, non-Abelian Bogomol'nyi-Prasad-Sommerfield (BPS) dyon pair at low energies. This system 
is amenable to the geodesic approximation which restricts the underlying SU(2) Yang-Mills-Higgs 
dynamics to an eight-dimensional phase space. We numerically calculate a representative set of 
long-time solutions to the corresponding Hamilton equations and analyze quasiperiodic and chaotic 
phase space regions by means of Poincare surfaces of section, high-resolution power spectra and Lya- 
punov exponents. Our results provide clear evidence for both quasiperiodic and chaotic behavior 
and characterize it quantitatively. Indications for intermittency are also discussed. 



I. INTRODUCTION 



o 
o 

The classical dynamics of non-Abelian gauge theories is known to be chaotic in a large part of its phase space Ul- 
^ ' By itself this is not unexpected since chaos is far more the rule than the exception in nonlinear dynamical systems. 
Perhaps more surprising, however, is the mounting evidence for this chaotic behavior, which is strictly speaking a 
classical phenomenon, to be of relevance for physical quantum gauge theories as well. 
I Part of the existing indications for the chaoticity of non-Abelian gauge theories stem from "homogeneous approxi- 
mations" which neglect all spatial variations of the fields. Although these drastic reductions of the dynamics access 
only a tiny and not generally physical fraction of the full phase space, the few remaining degrees of freedom proved 
sufficient to establish chaotic regimes first in SU(2) Yang-Mills theory |3| and later in Yang-Mills-Higgs theory Q and 
' Chern-Simons gauge theory Extensive lattice calculations of the time evolution under the full hyperbolic Yang- 
p\| i Mills equations subsequently showed that spatially inhomogeneous gauge fields not only evolve chaotically, too, but 
OO ' reveal even more complex and qualitatively new phenomena |6ll which are of direct physical relevance, for example, 

' in nonequilibrium processes |62l |. 
1^ I Between the above two computational extremes, i.e. either ad hoc truncations to only constant fields or the full 
solution of the classical field equations on the lattice, there exist physically interesting subsystems of the gauge 
dynamics whose spatially varying fields — typically soliton configurations — and time evolution can be studied without 
invoking uncontrolled approximations or requiring the solution of partial differential equations. In the following, we will 
Qh, focus on a prototy pica l such system, consisting of two electrically charged magnetic Bogomol'nyi-Prasad-Sommerfield 
^ ' (BPS) monopoles |l2l IT^ IT^. i.e. dyons, whose nontrivial spatial extension gives rise to crucial properties including 
the magnetic charge. Nevertheless, the low-energy time evolution of the dyon pair is accurately described by the 
geodesic approximation which involves just a few collective degrees of freedom, governed by ordinary differential 
equations. This geodesic dynamics is one of the best understood examples for classical (and quantum) interactions 
^ between extended solutions of physically interesting 3-1-1 dimensional field theories and can be formulated at a rare 
■ level of explicitness [Tsl H^. Beyond being fascinating in its own right (exhibiting e.g. the celebrated scattering 
angle tt/2 for head-on collisions), it therefore serves as a paradigm for the interactions among many other physically 
important solitons |6lj|. The purpose of the present paper is to examine regular and chaotic motion patterns of this 
system in a representative set of phase space regions. 

Besides appearing in the electroweak sector of standard model extensions, monopoles and their potentially chaotic 
interactions may have a crucial function in the context of quark confinement by the strong interactions. This becomes 
most explicit in the central role which BPS monopoles play in the confinement mechanism of A'^ = 2 supersymmetric 
Yang-Mills theory in 3-1-1 dimensions -By condensing in the vacuum, they realize the classic 't Hooft-Mandclstam 
dual superconductivity scenario |l8l | in which color magnetic charges get screened while color electric charges are 
confined by the dual Meifiner effect. Similar scenarios, in which the condensation of monopole-like objects plays a 
key role, are expected to unfold in more physical gauge theories as well [G^I- In 3+1 dimensional Yang-Mills theories, 
for example, there is lattice evidence for the condensation of Abelian-projected monopoles to generate the bulk of 
the string tension ■ (According to an interesting recent suggestion, the "active" monopoles might actually be 
BPS dyon constituents of caloron solutions with nontrivial holonomy |2ll|.) Hence monopoles ma y b e instrumental in 
resolving the most profound remaining mystery of QCD, i.e. the quark confinement mechanism |22l |. 

From a seemingly different but actually related perspective, quark confinement is expected to be linked to the chaotic 
behavior of classical chromodynamics as well. In fact, it has long been conjectured that the vacuum of non-Abelian 
gauge theories, when undergoing a transition from weakly to strongly coupled fields, also undergoes an order-disorder 
transition and that the strongly coupled QCD vacuum is populated by highly irregular color field configurations P]. 
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In the limit of a large number of colors, in particular, a vacuum made of random Yang-Mills fields has been shown 
to be a necessary and sufficient condition for quark confinement p^ . From the outset, one of the motivations for 
investigatingchaos in non-Abclian gauge theories was therefore to shed light on its potential role in the confinement 
mechanism Moreover, the instability of constant color-magnetic vacuum fields j2J| made it natural to conjecture 
that both gauge invariance and stability of the physical vacuum may be restored by disordering the color-magnetic 
background fields. Under the gluonic structures envisioned to carry the bulk of this disorder are random domains as 
well as populations of randomly distributed center vortices or monopoles. 

The last of these scenarios may be related to the subject of our investigation. Indeed, it is tempting to speculate that 
a disordered ensemble of monopoles (and antimonopoles) in a semiclassical vacuum may be generated by chaotic low- 
energy interactions among the monopoles. In the following we are going to investigate precisely this type of interaction 
in the simplest possible setting, i.e. between just two BPS monopoles. One may then hope that the quantitative 
understanding of its chaotic regime will lead to new insights into the disorder of interacting monopole ensembles as 
well. For sufficiently dilute systems, expansions in the monopole density and more sophisticated many-body techniques 
might even provide a starting point for the quantitative treatment of chaotic multimonopole ensembles. Alternatively, 
one could contemplate the technically challenging extension of the geodesic approximation to approximate BPS 
multimonopole-multiantimonopolc solutions which incorporate multimonopole interactions |65l |. 

Based on the above motivations, our main objective in the following will be to deepen the qualitative and quantitative 
understanding of quasiperiodicity and chaos in the geodesic motion of two BPS dyons. The pioneering numerical 
studies of this motion in Refs. I27L [281] already provided several indications for its nonintegrability. (Evidence 
for chaotic fluctuations around single monopoles in SU(2) Yang-Mills-Higgs theory exists as well although not 
too large, minimally spherically symmetric excitations remain regular ^^■) After a brief recollection of the classical 
two-dyon dynamics at low energies, we will first extend previous studies by examining Poincare sections for a set 
of numerically generated long-time trajectories of the dyon pair. In the subsequent sections we break new ground 
by invoking high-resolution power spectra and maximal Lyapunov exponents to analyze the dyon orbits further. 
This analysis will go beyond the mere identification of standard motion patterns and provide the first quantitative 
characterizations of quasiperiodic and chaotic two-dyon trajectories. 



II. CLASSICAL LOW-ENERGY DYNAMICS OF THE TWO-DYON SYSTEM 



In order to set the stage for our investigation, we first recapitulate some pertinent aspects of the dynamics of a 
classical BPS dyon pair at low energies and discuss its known constants of the motion. Readers interested in more 
detail are referred to the lucid discussion by Gibbons and Manton (15>] . 

We consider Yang-Mills-Higgs (YMH) theory in 3+1 dimensions with gauge group SU(2) and the Higgs field in 
the adjoint representation, i.e. the Georgi-Glashow model. The classical field equations admit topological soliton 
solutions which are magnetic monopoles with integer magnetic charge k [3l| . In the following we will be interested in 
these monopole solutions in the BPS limit of vanishing Higgs potential |l2l[l3lll3|. which solve the more restrictive 
Bogomol'nyi equation 

= ieyfcF/, = ± (j-^a,; + .96""=^^) (1) 

(where F^^^ is the field strength tensor of the gauge field ^4° and is the (adjoint) Higgs field). The solutions of 
Eq. are the absolute minima of the YMH energy in their topological charge sector and form a submanifold Mk 
in the space of gauge-inequivalent finite-energy fields. Static multimonopole solutions (with |fc| > 1) of Eq. ^ are 
possible because the repulsive magnetic forces between the individual monopoles are counterbalanced by the attractive 
forces which the massless Higgs field mediates |3^ . 

Although the underlying YMH dynamics is not directly embedded in the standard model, it appears naturally in 
grand-unified theories. Monopole solutions similar to the BPS prototype might therefore be physical. Their mass 
would probably be large enough to explain why they have so far escaped discovery in earthbound laboratories. 

The monopole solutions come in families whose members are characterized by continuous collective coordinates or 
"moduli" x" . (For the one-monopolc solution, for example, these arc the three position coordinates of the center and 
an overall phase angle.) Hence the moduli space spanned by the collective coordinates is just the (generally curved) 
manifold Mk- Its metric gap^x) is induced by the metric on the more comprehensive space of all finite-energy field 
configurations which the kinetic terms in the YMH Lagrangian define. 

The time evolution of a nonstatic and therefore in general electrically charged BPS monopole system is governed by 
the hyperbolic partial differential YMH equations whose quantitative analysis and solution poses rather formidable 
problems. Reassuringly, however, the low-ene rgy dynamics of BPS dyons reduces to a much more tractable problem 
as long as Manton's geodesic approximation j33l | can be invoked. The latter rests on the observation that energy 
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conservation forces the low-energy motion of the fc-dyon system (with sufficiently small initial velocities of the collective 
coordinates, tangent to M^) to remain close to one of the static BPS solutions on Afj. at all times. This is simply 
because contains all absolute YMH energy minima with magnetic charge fc, so that moving out of it would cost 
both kinetic and potential energy. Hence the low-energy motion of a fc-dyon system approximately corresponds to 
low-energy motion of an associated point on the moduli space M^. Since the energy of all static fc-monopole solutions 
is degenerate, i.e. at the same (minimal) potential, furthermore, the low-energy motion is approximately determined 
by the kinetic energy alone. Hence it corresponds to geodesic motion on the moduli space Af^, which is governed by 
the purely kinetic Lagrangian 

Lgeod = ^gap{x) (2) 

where m is the reduced mass of the dyons. Physically, this just means that at small velocities (compared to the 
velocity of light) internal excitations (vibrations) and deexcitations (radiation) can be neglected, i.e. the dyons adapt 
adiabatically to their interactions by deforming reversibly and scattering elastically. (Corrections to the geodesic 
approximation were analyzed in an effective field theory framework in Ref. |34j.') 

In the following, we will focus on the two-dyon system. Because of the product structure of the moduli space, its 
center of mass momentum and an overall phase (whose time dependence is associated with the total electric charge) 
are individually conserved, and their metric is flat. Hence those degrees of freedom decouple from the internal motion 
and can be separated out. The remaining dynamics simplifies to the geodesic motion in the four-dimensional internal 
part M^'' of the moduli space and can be studied independently. A physically intuitive coordinate system on M^^ 
consists of the Euler angles i9, and '0, which determine the orientation of the two-dyon system [b^ . and the distance 
variable g which measures the separation between the two dyon centers (at large g). 

Although the metric on M^^ is induced by the SU(2) YMH dynamics, its direct derivation from the kinetic terms 
of the YMH Lagrangian seems out of reach. Instead, it has been constructed on the basis of ingenious symmetry 
arguments by Atiyah and Hitchin (AH) ^| , as summarized in Appendix ^ Specialization of Eq. ||2J) to the AH 
metric then determines the geodesic dynamics of the dyon pair explicitly. The resulting Lagrangian 

Lah = \ (? + a\g) ujI + h\g) c,^ + c\g) c.^] (3) 

turns out to be nonlinear since the AH metric is curved. The functions a, 6, c and / of the separation g are given in 
Appendix IXI and the Wi('i?, v?, ^Z") are angular velocities of the monopole (or dyon) pair around the axes of the body- fixed 
frame, 

LOx = sin ^ + ip sin -d cos ip , (4) 
LOy = d cos t}} + sin d sin tp , (5) 
= -0 + i^cos?? . (6) 

The form of the Lagrangian (jS)) is analogous to that of a nonrigid body with distance-dependent "moments of inertia" 
a? and c? around the body-fixed axes. Following Gibbons and Manton, we define the radial coordinate g by 
choosing / = —b/g which leads to convenient expressions for a, b and c |15|. In Appendix^ the four Euler-Lagrange 
equations of the geodesic motion are derived by variation of Eq. Ipjjl. 

For the question of integrability versus chaos, the number of integrals of the motion plays a decisive role. In fact, 
a motion is (Liouville) integrable only if the number of independent conserved quantities (whose mutual Poisson 
brackets vanish) at least matches the number of degrees of freedom. For the geodesic dynamics (jSJ , three constants 
of the motion are known explicitly (cf. Appendix IB 2|l . namely, the total angular momentum 

= - 2p^p^, cot z9 esc z9 + (p^ -f- p^) csc^-d (7) 



(where the Pa are generalized momenta canonically conjugate to the coordinates a, as defined in Eq. IIB3|) of Appendix 
r5)l. the energy ()B17|) which, for geodesic motion, equals the Lagrangian since the potential on Afj^"* vanishes. 



EAH^^\^+a'ul + b'c.l+cW\ (8) 




and, finally, the generalized momentum p^p conjugate to the coordinate (p which is cyclic, i.e. does not appear explicitly 
in Eq. @. 
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At least one additional, independent constant of the motion is therefore required for the two-dyon motion to become 
integrable. Such a fourth conserved quantity indeed exists (to arbitrarily good approximation) in at least one region 
of phase space, namely, where the two dyons remain infinitely separated. Since they cannot exchange electric charge 
then j6^, the time evolution of each of their phases or, equivalently, their individual electric charge is conserved. 
Hence the motion of two far separated dyons must be integrable, as intuitively expected, and cannot exhibit chaos. 
(At asymptotic distances the AH metric reduces to the Euclidean Taub-NUT metric whose geodesic motion is indeed 
known to be integrable [TsIIb^.) 

This situation changes, however, if the two dyons begin to approach each other. Only their total charge, but not 
their charge ratio, remains conserved since the dyons are then able to exchange charge through the Higgs field. (Even 
when starting asymptotically with two uncharged monopoles, one will therefore generally end up with two dyons of 
opposite but nonvanishing charge.) If the integral of the motion associated with the relative phase ceases to exist, 
chaotic motion becomes possible |68j . It is this region of the phase space in which we will be particularly interested 
below. 

The above observations raise the question for which orbits the violation of individual charge conservation can 
typically be considered as a small perturbation away from the integrable Taub-NUT limit. In our context, this 
question is of relevance because for such orbits the transition from integrable to chaotic motion will be delayed and 
potentially obscured by the implications of KAM theory. Indeed, the KAM theorem [s^l states that almost all invariant 
tori of the unperturbed Taub-NUT motion will remain intact during sufficiently small deviations from the g oo 
limit 6Qi|. Under weakly nonintegrable perturbations the motion should therefore stay quasiperiodic for almost all 
initial conditions. Chaotic motion will then be restricted to the small set of trajectories which lie on the descendants 
of invariant tori with commensurate frequencies. 

Although we have argued above that it is unlikely for the geodesic AH motion to be integrable outside of the 
asymptotic region, the existence of additional but so far undiscovered constants of the motion cannot be excluded a 
priori since no general method for finding all conserved quantities of a nonlinear dynamical system is available. Hence 
it was an important step by Temple-Raston and Alexander to gather the first numerical evidence for the existence of 
chaotic regions in the two-dyon phase space '2^, ^ 271 [2p . In the next section, we will elaborate on part of these results 
by extending the Poincare section analysis of Ref. [23 and by looking for orbits which bear the insignia of chaos. 

III. POINCARE SECTIONS 

We begin our search for chaotic regions in the phase space of the dyon pair by analyzing the Poincare sections of 
several typical trajectories. A Poincare section draws a selective portrait of a given phase space orbit (over a finite 
time interval) which achieves an enormous data reduction by mapping the whole trajectory into a discrete set of 
points. This set transparently exhibits characteristic global aspects of the orbit and allows, in particular, a direct 
visual distinction between orbits which arise from integrable (or weakly nonintegrable, in the KAM sense) and fully 
nonintegrable dynamics. 

Poincare section analyses are most powerful for conservative systems with two degrees of freedom where a trans- 
parent graphical interpretation becomes possible. At first glance one might therefore doubt their utility in the eight- 
dimensional phase space of the two-dyon problem. However, in Ref. it was recognized that a particular canonical 
transformation turns a second constant of the dyon pair motion (besides the generalized momentum p^p), namely, the 
total angular momentum squared M^, into a canonical momentum. As a consequence, the transformed Hamiltonian 
becomes independent of the associated canonical variables while and act as fixed "external" parameters. This 
Hamiltonian actually defines a reduced, four-dimensional (and still symplectic) phase space in which p^ and 
are automatically conserved at every point. Energy conservation further constrains all orbits to a three-dimensional 
hypersurface inside this reduced phase space and thus makes a graphical Poincare section analysis feasible. 

The Poincare section of any given orbit is the set of its intersection points, in a fixed direction, with a suitable 
two-dimensional plane in phase space. In our case, a useful choice is the {g,Pg) plane located at the position 
Ml = which we will adopt throughout. Hence we define the elements of an orbit's Poincare section as those points 
in which it pierces through this plane while Mi changes from positive to negative sign. 

One can easily convince oneself that integrable and chaotic motions generate qualitatively different Poincare sections. 
Indeed, the existence of an additional, fourth constant of the motion C4 (as provided e.g. by p^ for far separated 
dyons, cf. Appendix IB 2\i would restrict all orbits to a generally two-dimensional submanifold in phase space, namely, 
the intersection of the two three-dimensional hypersurfaces on which either the energy or C4 have fixed values. 
The dyon-pair motion then becomes (Liouville) integrable, the submanifold becomes an invariant torus and the 
orbit's intersection points with the {g,Pg) plane sweep out a (maximally) one-dimensional curve. If the system is 
nonintegrable, on the other hand, no fourth constant of the motion exists and the orbit is only bound to the three- 
dimensional constant-energy surface (for strong enough nonintegrability in the KAM sense). Its intersection points 
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with the {q,Pq) plane can therefore be more broadly distributed and eventually fill a two-dimensional area, which 
provides a clear signature for a chaotic dynamical regime |70| . 

In order to prepare for the Poincare section analysis, we have generated twelve representative orbits by numerically 
integrating the equations of motion IIB12|l - IIB15|l under suitable initial conditions |7ll |. The latter are provided by 
specifying initial values for the four coordinates and their time derivatives. In order to facilitate the comparison with 
the results of Ref. |27j . we adopt the following five initial conditions, and the form of a sixth, at the initial time t^: 

t}=^, g = ^ = and ip^^. (9) 

As a consequence of Eqs. (jB3ll . the corresponding initial values of the conjugate momenta are 

Pe = 0) Pi) ^M2^ 0, = Ml = /i2 and = M3 = c^-ip, (10) 

which implies 

= MI+MI+MI=pI+pI+ pI = hl + cV- (11) 

Here /12 is the parameter in the reduced Hamiltonian which fixes the value of the conserved momentum p^p. Hence only 
the initial values for ij: and g remain to be determined. If /12 is given, Eq. Hll|) implies that the initial condition for 
can be specified by fixing the value of another conserved quantity, namely, the total angular momentum squared . 
Finally, instead of prescribing the initial condition for g directly, it is more convenient to use the reparametrization 

^? = 2A-(^sin0 (12) 

{K is an elliptic integral, cf. Eq. IjAlSp ). which maps infinite dyon separation {g 00) into /3 = tt, and to specify the 
initial value of p. Under the above conditions (O, each orbit is therefore uniquely characterized by the values of /3 
and two integrals of the motion, /12 and M^. For the interpretation of the results below it will be helpful to keep in 
mind that (at fixed /12) the initial relative electric charge (to) of the two dyons, and therefore their initial Coulomb 
interaction, grows with (cf. Eq. Ullfl '). 

A fourth-order Runge-Kutta algorithm |39| | was used for the numerical integration of the equations of motion. The 
Eqs. (|A14() are solved along the way to obtain the values of the AH functions a, b and c at each g. The length of an 
integration step was typically 10"'' and the longest runs contained 2^^ steps. The accuracy of the generated orbits 
was monitored after each time step by calculating the three conserved quantities (i.e. E, 3Vp and p^). The maximal 
deviations from their fixed values were of the order 10~^^. The Poincare section of a given orbit is then constructed 
from those of its points at which Mi vanishes and Mi < |72j . Since the numerical integration routine generates the 
orbits only at discrete time steps, the exact intersection points were determined by polynomial (Lagrange) interpolation 
between the two adjacent points on the orbit with opposite signs of Mi. 

We will now discuss the Poincare sections for twelve orbits which all share the same /12 value, /12 = 37.596. We 
arrange them into three plots, each containing the sections of four orbits with a common value of M^ and initial 
values P = (3.11,3.12,3.13,3.14) of the radial coordinate. The corresponding initial values for {g,Pg) are indicated 
by a small circle although Mi = ft.2 7^ implies that they do not lie in the Mi = plane. Each of the twelve selected 
orbits describes charged monopole, i.e. dyon interactions since M^ > = 1413.5 and Eq. Hll|) imply that the initial 
value of their momenta = c^?/', associated with the relative electric charge, is nonzero. The above initial conditions 
were chosen to yield a rather representative set of Poincare sections similar to those considered in Ref. , in order 
to both qualitatively confirm the results obtained there and to extend them into neighboring phase space regions. 

The first plot, Fig.^ contains the four Poincare sections with the smallest value M^ = 1906.71, i.e. with the weakest 
initial Coulomb attraction between the dyons. This explains why each of the orbits covers a rather large range of g 
values while the variation of the radial velocities, i.e. the Pg range, is relatively moderate. In three of the sections, 
furthermore, the two dyons never come closer than their initial distance, and in all four their separation g ^ 11 seems 
to stay well inside the asymptotic region where the relative charge of the dyons becomes (almost) time-independent 
and the motion (approximately or KAM-) integrable. This is confirmed by the Poincare sections whose points indeed 
trace a one-dimensional closed curve which, within plot resolution, appears mostly continuous. Poincare sections of 
this type are generated by superpositions of periodic motions with incommensurate frequencies, i.e. by quasiperiodic 
orbits. As discussed above, quasiperiodic behavior indicates that the dynamics is either integrable or at most weakly 
nonintegrable in the KAM sense. 

In the next plot. Fig. |21 we display the analogous Poincare sections with a larger value of M^ ~ 2152.95 and 
therefore with a stronger initial Coulomb attraction. Comparison with the sections of Fig.^shows that the variations 
in dyon distance are now smaller (tighter orbits) while their relative momenta vary more strongly over each of the 
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orbits. The increased attraction also brings the dyons closer together, their maximal separation now being the initial 
one for all four orbits. Nevertheless, their minimal separation Qmin ~ 8 seems to stay large enough for the motions to 
remain approximately integrable since all four Poincare sections still form one-dimensional, closed curves. 

We therefore increase further, to = 2237.70, and plot the corresponding Poincare sections in Fig.O Clearly, 
the character of the outermost section differs qualitatively from all those encountered previously. Instead of remaining 
constrained to a curve, it visibly spreads out into a two-dimensional area of the {g,Pg) plane. Hence it corresponds 
to an aperiodic orbit and strongly suggests that the dynamics has become chaotic. This interpretation is consistent 
with the fact that the two dyons approach each other the most in this orbit. Their minimal distance gmin ~ 27r, which 
is about twice the distance of the bolt singularity at p = tt, apparently suffices for electric charge exchange by the 
Higgs field to become efficient and hence for to cease being even an approximate integral of the motion |73 | . This 
confirms earlier indications for the chaoticity of the geodesic two-dyon dynamics |26l | which were supported by the 
analysis of similar Poincare sections [23|, Julia sets ^3 ^^id escape plots of two-monopole scattering trajectories [28l |. 
An additional chaotic orbit, with a still larger value of M^, will be generated and analyzed in the following sections. 



IV. POWER SPECTRA 



Numerical studies of a single chaos indicator are generally not sufficient to establish chaotic behavior with certainty. 
Inevitable numerical roundoff errors as well as specific features of a system under consideration (e.g. particular regions 
of instability) make it often desirable to probe the character of the motion by several complementary techniques. 
Typical pitfalls, like the premature misinterpretation of a seemingly irregular but perfectly integrable quasiperiodic 
motion pattern as chaos, can thereby be avoided. 

Power spectral analysis has proven particularly useful for the distinction between (quasi-) periodic and chaotic 
time evolution • This is because the power spectra of ordered motion (either periodic or quasiperiodic) consist of 
sharp resonance lines which appear at simple harmonics of the base frequency in the periodic case and at any linear 
combination of all integer multiples of the base frequencies in the quasiperiodic case. Aperiodic systems, in contrast, 
are generally chaotic and have continuous and noisy power spectra. 

In order to sharpen and extend the interpretation of the Poincare section analysis of Sec. IIIII we have submitted 
the underlying orbit data to a spectral analysis. These data are four-dimensional, discrete time series of specific 
orbit solutions on Mj"' whose coordinates {x"{tk)} = {'d,(p,ip: 9} &re recorded numerically at equally spaced times 
tk := fcAt. Based on these solutions, one can calculate the time evolution of any dynamical variable f{{x'^{tk)}) of 
interest as well as its discrete Fourier transform. We found it useful to choose the generalized momentum Pg(tk) ~ 
f'^{g(tk)) g{tk) conjugate to g as the dynamical variable. It has the power spectrum 

2 

(13) 

which is a function of the mode frequencies h'. (The Wiener-Khinchin theorem asserts that the alternative definition 
of P{iy) as the Fourier transform of the time scries' autocorrelation function is equivalent if the correlations decay 
sufficiently fast.) In order to assess a potential bias due to the choice of pg as the dynamical variable, we have also 
calculated the power spectra associated with the time evolution of g and obtained essentially analogous results. 

We numerically compute the power spectra by means of the Fast Fourier Transform algorithm. Appropriate 
windowing is used to suppress artificial oscillations due to the discrete data set and the finite time interval (see 
Ref. a-iid references therein). Since the impact of these artifacts decreases with the length of the time interval 
over which data for Pg{tk) are available, we sample the orbits after each of 2^^ time steps of length At = 10~* which 
ensures sufficient accuracy for our purposes. We have calculated power spectra for the three orbits associated with 
the outermost Poincare sections in each of the Figs. specified by the initial conditions given in Sec. IIIII and for 
an additional orbit to be discussed below. 

The Pg{tk) time series based on the orbit associated with the outermost Poincare section of Fig. ^is depicted in 
Fig. EJi over 4 x 10^ time steps. The corresponding power spectrum (obtained from all data on pg, i.e. from the full 
orbit with 2^^ time steps over a total time of 3.35 x 10'^) is shown in a relatively small but representative frequency 
interval in Fig. On the basis of its Poincare section, the underlying orbit was interpreted as quasiperiodic in 
Sec. IIIII The power spectrum confirms this interpretation but contains much more quantitative information. Indeed, 
while the underlying time series in Fig. 0^ is in many ways indistinguishable from a periodic one, its power spectrum 
exhibits sharp peaks |74| of varyin g st rength which are located at odd integer linear combinations of two rationally 
independent (or incommensurate) [75j base frequencies i^i ~ 0.0653 and 1^2 — 0.0978. Hence this orbit is two-mode 
quasiperiodic [t^. 
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We now turn to the orbit whose Poincare section is the outermost in Fig. |21 Again, its time series (cf. Fig. ^Jd) 
is visually difficult to distinguish from periodic motion while its power spectrum in Fig. 2f reveals a quasiperiodic 
motion with the two incommensurate base frequencies I'l ~ 0.1234 and 1^2 — 0.1848. In contrast to the previous 
power spectrum, however, the peaks at frequencies corresponding to more complex linear combinations of i^i and 
1^2 are practically undetectable here. This seems to be a consequence of the larger initial coupling between the two 
dyons and is a rather frequent occurrence among sufficiently strongly interacting quasiperiodic systems. Although 
experimentally confirmed in many situations, a rigorous theoretical explanation for this type of behavior appears still 
to be missing |4^ . 

As foreshadowed by the results of the Poincare section analysis, a strikingly different power spectrum is obtained 
from the orbit associated with the outermost section of Fig. |21 Indeed, this is the orbit whose Poincare section 
indicated the chaoticity of the geodesic dynamics. Although the time dependence of its Pg in Fig. 01: remains, apart 
from small and irregular modulations, similar to the previous cases, the corresponding power spectrum in Fig. 01; has a 
qualitatively different character. Instead of sharp and isolated peaks, it now contains equally spaced, broadened peaks 
on top of a smooth and (within our frequency resolution) continuous background power distribution , covering all 
recorded frequencies and signalling the onset of aperiodicity . 

The obvious departure from quasiperiodic behavior strengthens the evidence for the chaoticity of this orbit. Our 
discussion in Sec. ^1 indicates, furthermore, that even the quasiperiodic behavior of the previous spectra, especially 
that of Fig.0f, may have been generated by orbits which lie outside of the strictly integrablc phase space region corres- 
ponding to asymptotically far separated dyons. Indeed, a discrete power spectrum still results if the nonintegrability 
is caused by weak perturbations and if the orbit remains on a KAM torus, i.e. continues to behave quasiperiodically 
|79j . Our results imply that such a "delayed" onset of chaos might persist down to minimal distances gmin ^ 8 between 
the two dyons. 

In any case, the above findings suggest that the faster and closer the two dyons approach, the more unpredictable, 
i.e. chaotic their motion becomes. In order to test and extend this conclusion, we have calculated a further orbit with 
a still larger total angular momentum squared (and consequently stronger initial Coulomb interaction), = 2359.46, 
and a somewhat smaller initial distance f3 = 3.13 between the dyons than in the previous chaotic orbit. Its power 
spectrum, displayed in Fig. ^Ji, has a dense background and confirms the expectation that the new orbit is aperiodic 
as well. Comparison with the power spectrum of the first chaotic orbit shows, furthermore, that the background 
in Fig. 2Ji is less noisy and that more sharp frequencies remain clearly discernible than in Fig. The faster time 
dependence of the associated Pg(t) in Fig.^, on the other hand, exhibits more pronounced amplitude modulations. 
In the following section we will continue the analysis of the two chaotic orbits by calculating their maximal Lyapunov 
exponents. 

One might perhaps wonder why we have encountered no quasiperiodic behavior with more than two fundamental 
modes. This is not particularly surprising, however, since two-mode quasiperiodicity is the rule rather than the 
exception among sufficiently strongly coupled, nonlinear dynamical systems. Indeed, the nonlinear couplings between 
an increasing number of modes tend to replace quasiperiodicity by chaos . It is quite plausible that this happens 
in the two-dyon system as well and explains the predominance of quasiperiodic orbits with the minimal number of 
two base frequencies. 

In summary, our spectral analysis has uncovered that those parts of the two-dyon phase space which remain close 
enough to the asymptotic g 00 region are characterized by two-frequency quasiperiodic motion. In addition, we 
have determined the base frequencies of two typical quasiperiodic orbits quantitatively. Perhaps most importantly, the 
power spectra have also accomplished their main task to separate quasiperiodic from irregular behavior. In particular, 
they strongly support the identification of two orbits, the last one from Sec. IIIII and an additional one with a still 
larger initial Coulomb force between the dyons, as chaotic. Our findings therefore substantially increase previous 
evidence that, apart from the asymptotic g 00 region, the relative low-energy motion of two BPS dyons admits 
only three independent conserved quantities and turns out to be genuinely nonintegrable. 

While the distinction between (quasi-) periodic and chaotic motion is a particular strength of the spectral analysis, 
it does relatively little to further characterize chaotic behavior [s^l- Hence we consider it useful to subject our orbits 
to yet another classic analysis tool from the arsenal of chaos indicators, by calculating their characteristic or Lyapunov 
exponents. These exponents complement our previous analyses particularly well since they are specifically designed 
to quantify the chaoticity of irregular motion patterns. 

V. LYAPUNOV EXPONENTS 

While the high-resolution spectral analysis of the last section has very clearly distinguished quasiperiodic from 
irregular orbits and determined the type of quasiperiodic behavior and its base frequencies quantitatively, it is much 
less specific about the properties of the aperiodic motion. It cannot, for example, unambiguously distinguish chaotic 
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from (quasi-) random behavior (potentially due to roundoff errors). 

Although previous work and our results of the last two sections have produced strong evidence for the chaoticity of 
two-dyon orbits in specific regions of the (finite g) phase space, it therefore remains desirable to carry the investigation 
a step further and to obtain a quantitative characterization of the chaotic behavior. The most fundamental such 
characterization is provided by the values of positive Lyapunov exponents |47j . After a brief and informal introduction 
to the underlying concepts, we will therefore evaluate the largest of these characteristic exponents for selected dyon 
pair orbits. 

The primary function of Lyapunov exponents is to quantify the logarithmic rate of convergence or divergence 
between two orbits which started at an initial time to at neighboring positions Xg and xq + Sxq (with ||Jxo|| <C 1 
where ||-|j denotes the norm with respect to a Riemannian metric). Their relevance is obvious: exponentially divergent 
orbits, corresponding to positive Lyapunov exponents and signalling an exponential sensitivity of the time evolution 
to the initial conditions, are the prototypical signature of chaos. More specifically, under dynamical evolution for a 
sufficiently long time t the deviation between the position of both orbits becomes Jx^ and the maximal Lyapunov 
exponent Lmax is defined by relating the deviation norms as \\S^t\\ ~ exp(L,„axi) ||'5xo|i, i.e. 

Uax= lim lim jln^^. (14) 

t-»oo ||5xo|H0 t ||c)Xo|| 

For a detailed and rigorous treatment of Lyapunov exponents we refer to original work by Oseledec ^3 ^^i^ Ruelle 
[4^ . as well as to the review |4^ . 

Although the definition 1)14(1 of Lmax is conceptually transparent, it does not lend itself to direct numerical imple- 
mentation since in chaotic systems any initial deviation ||5xo|j, no matter how small, will eventually evolve into a 
number ||^Xi|| which far exceeds the (floating point) representation capabilities of computers. Hence more indirect 
numerical approaches are called for. and several different ones have been developed over the last decades jSOj . In 
the following, we will adopt the Jacobian method which integrates the equation for the time evolution of the devia- 
tion between two initially neighboring orbits, linearizing it anew at each time along the orbits and therefore having 
to manipulate only the relatively small deviations produced by one time step. A generalization of this method to 
the calculation of the whole non-negative Lyapunov spectrum was developed in Refs. |5ll l53 |. Our discussion and 
implementation follows Ref. [s^ . 

The Jacobian method is most transparently formulated by rewriting our system l|B12p - l|B15p of four second-order 
differential equations into a system of eight autonomous first-order equations 

i(i) = F(x(i)), (15) 

where x comprises the four coordinates and their time derivatives, i.e. x and F are eight-dimensional column vectors. 
Under the initial conditions x{to = 0) = Xq, the unique solution of Eq. (|15|l is the orbit x(t) = 4>(i,xo) with xq = 
*(0,xo). We now linearize the above system in the small deviation dx{t) at any point x(t) of the trajectory by 
substituting x(f) = x{t) + Sx{t) into Eq. (|15|l and neglecting terms of second-order in 5x{t). The result is a linear 
system of first-order equations 

Sxit) ^ 3{x{t))-Sx{t) (16) 

for the deviation 6x which contains the Jacobi matrix J(x(t)) ~ [(?F(x)/(?x]|x=x(t)- In terms of the tangent vector 
^{t) = limj|5x(o)||-+o <^x(0/ll'^^(0)ll the orbit at x{t), Eq. (|16() turns into the so-called variational system for ^{t) 
under the initial data xq, 

m = Ji'i>{t,xo))-m- (17) 

It propagates small variations tangent to the orbit at time tg to small variations tangent to the orbit at time t. In 
terms of the ^, Eq. H14|l for L„iax becomes 

Lmax(xo,^o)- lim X(t) where ;^(t) := i^^MM. (18) 

t^oo t 

Hence, one can calculate Lmax by numerically integrating the variational Eqs. (|17|) up to sufhciently large times. In 
order to supply the necessary input, i.e. the Jacobian J($(t,xo)) at each time step, one has to solve the equations 
of motion (|15|l for the orbit x{t) = ^{t,xo) in parallel. A biased choice for the initial tangent vector is avoided by 
selecting its orientation randomly. The long-term time evolution of the orbits and ^(i) will then be dominated by the 
largest Lyapunov exponent |8lj |. 
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The numerical evaluation of Eq. (|18|l requires some additional precaution, however, since for chaotic orbits the norm 
will grow large enough to generate floating point overflows on computers. To circumvent these, we directly 
calculate the value of the required logarithmic norm after n time steps (i„ = nAt) as the sum 

ln!lC(t„)!l=X:ini|^ (19) 

(with ||^(to)|| = ll^oll = 1) '^^ logarithmic length increments after each time step [53 |. The incremental norm changes 
||^(ii)||/||4(ti_i)|| remain small enough to be representable by floating point numbers and are obtained by renormalizing 
^ to unit length after each evaluation. (Of course, the renormalized ^ remain solutions of the linear variational system 

We now calculate, by means of the Jacobian technique, the maximal Lyapunov exponents of the four dyon-pair 
orbits whose power spectra were obtained in the last section and plotted in Figs.Qle-h. The Poincare sections of the 
first three of these orbits are the outermost in Figs. EE] The integration range consists of n = 2 x 10^ time steps 
which corresponds to approximately 650 periods for the outermost quasiperiodic orbit of Fig. ^ In order to monitor 
the time evolution of x during its approach to the limiting value L^ax according to Eq. (|18|l . we plot x{t) for each of 
the four orbits in Fig. [3 Obviously, there is a striking qualitative difference between the plots for the two orbits which 
we had identified as (two-frequency) quasiperiodic in Sec. IIVI (with power spectra in Figs. 2)3 and^f) and the irregular 
ones whose power spectra are shown in Figs. 0^ and For the quasiperiodic orbits one infers within numerical 
uncertainties that L^ax = IhTit^ooX(i) ~ while x{t) approaches the finite and positive values 

Lmax,l 0.02, Uax,2 ^ 0.008 (20) 

for the two aperiodic orbits. Around these orbits the motion is therefore exponentially sensitive to small variations of 
the initial conditions in at least one phase-space direction. This is the prototypical hallmark of chaos. We have thus 
achieved two of our main objectives, namely, the unequivocal confirmation of the chaoticity of the two-dyon system 
and a (semi-) quantitative determination of its primary characteristic scales. 

The maximal Lyapunov exponent L^ax,! of the first chaotic orbit (with the power spectrum in Fig. ^) is more 
than twice as large as that of the second chaotic orbit (with the power spectrum in Fig. 0J.i), L,nax,2, although the 
differences in the initial conditions (M^ = 2237.70 and /3 = 3.14 for the first and = 2359.46 and /? = 3.13 
for the second orbit) seem comparatively small. Closer inspection of these two orbits reveals that the minimal and 
maximal dyon distances in the first one are Pmin,i = 6.99 and gmax.i — 17.98 while the corresponding momenta pg 
vary inbetween Pg,imin/max = ±2.09. For the second chaotic orbit one has a smaller variation between minimal 
and maximal g values, f3min,2 = 6.73 and £imax,2 = 14.23, and a somewhat smaller range of momenta, bounded by 
Pe,2min/max = ±1-94. Although the initial Coulomb attraction between the dyons is stronger in the second orbit (since 
its total angular momentum squared is 5% larger) and the dyons therefore come closer to each other and depart 
farther from the integrable asymptotic domain, its chaoticity — as measured by the maximal Lyapunov exponent — is 
still considerably smaller. This might be a consequence of the generally smaller momenta of the second orbit. 

The behavior of x also reveals a conspicuous qualitative difference between the two chaotic orbits. While the x{t) of 
the first one (dot-dashed curve with filled triangles in Fig. (Sj) clearly stays above those of the quasiperiodic orbits for 
all t, the x{t) of the second one (dashed curve with filled squares in Fig.O follows its quasiperiodic counterparts for a 
long time rather closely and then suddenly rises in a "burstlike" onset of chaos. This intriguing behavior may be a first 
vestige of intermittency in the geodesic dyon-pair motion. It also seems to explain the more regular and approximately 
quasiperiodic power spectrum of the second orbit [Fig.^]. Stronger evidence for the potentially intermittent behavior 
could be established by identifying complete intermittency intervals of x(t) but would require much longer orbits. 
The power spectra of Sec. IIVI provide further testing grounds for intermittency which could be exploited, e.g., by 
subjecting them to a moment analysis or by searching for a power-law behavior in their low-frequency tails (i.e. l/f 
noise) |5q| . 

Although the large-time limit of x(i) (i.e. the maximal Lyapunov exponent) vanishes for quasiperiodic orbits, its 
characteristic time dependence may contain useful information as well. In order to examine the behavior of x{t) for 
quasiperiodic motion more closely, we select a new sample of three quasiperiodic orbits, namely, those whose Poincare 
sections are the mnermost in each of the Figs. and display the corresponding x(t) in Fig. |B1 (The full (dotted, 
dash-dotted) line corresponds to the orbit with the innermost Poincare section in Fig. ^ EJ-) At sufficiently large 
t all curves seem to approach straight lines, which indicates a power-law behavior 

X(t)^at-'' {a,b>0). (21) 

This type of scaling behavior is frequently encountered in quasiperiodic systems and supports the visual impression 
that all curves have indeed the expected L,„ax = liint^ooX(0 = 0- 
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We close this section by recalling that (positive) maximal Lyapunov exponents contain crucial information about 
the physical behavior even of quantum systems. A notable example is their partial characterization of nonequilibrium 
processes in semiclassical systems (e.g. at high temperature) where they typically set the scale of relaxation times 
and thcrmalization rates [ll|. One might expect that our Lyapunov exponents play a similar role in determining, e.g., 
the equilibration rate of a nonequilibrium dyon system. 

VI. SUMMARY AND CONCLUSIONS 

We have analyzed several representative motion patterns of two interacting BPS dyons in the geodesic approxi- 
mation. The main emphasis was put on discerning regular and chaotic orbits and on characterizing them both 
qualitatively and quantitatively by means of suitable chaos indicators. 

Our study is based on a sample of thirteen long-time phase space trajectories for which four-dimensional time series 
were generated by numerically integrating the equations of motion with high accuracy over typically 2^^ time steps. 
The initial data sets were chosen to cover a representative range of motion patterns and to explore the low-energy dyon 
interactions at different strengths. Hence the orbit set includes sequences of trajectories whose decreasing minimal 
dyon separations interpolate between asymptotic dyon distances, where charge exchange becomes ineffective and the 
geodesic dynamics integrable, and relatively small minimal separations for which the interactions are expected to 
become nonintegrable. 

A second motive for our initial-data selection was to include orbits in the vicinity of those for which Poincare 
sections were already available. This allows for a direct comparison with previous results and served as a useful 
benchmark and starting point for our work. We constructed Poincare sections in the radial coordinate-momentum 
plane for twelve orbits. They qualitatively confirm the earlier results and extend them to trajectories in neighboring 
phase space regions. Moreover, they contain useful graphical information on the shape of the constant-energy surfaces 
and on the qualitative behavior of the two-dyon system as a function of the initial conditions. The dimensionality of 
the Poincare sections, in particular, provides clear indications for the underlying trajectories to be either quasiperiodic 
or chaotic. 

In order to complement the results of the Poincare section analysis, we have additionally calculated high-resolution 
power spectra of selected orbits. The spectral analysis of the momentum conjugate to the dyon separation has 
provided particularly clean distinctions between quasiperiodic and aperiodic orbits which strengthen and extend the 
interpretation of the Poincare sections. In addition, the power spectra produced the first quantitative characterization 
of quasiperiodic dyon-pair orbits by establishing the number of their fundamental modes (two), determining their 
frequencies and yielding the strength distribution over the various harmonics. The emergence of just the minimal, i.e. 
two-mode quasiperiodicity is rather widespread among sufficiently strongly coupled, nonlinear dynamical systems. 
The common expectation that nonlinear couplings between more than two fundamental modes increasingly turn 
quasiperiodicity into chaos might therefore apply to the two-dyon system as well and explain why we have only found 
two-modc-quasiperiodic and chaotic trajectories. 

In contrast to their almost complete characterization of quasiperiodic motion patterns, power spectra do extract 
relatively little pertinent and quantitative information from irregular orbits. Hence we have additionally calculated 
the primary characteristic scales of chaotic motion, i.e. the maximal Lyapunov exponents, for a suitable subset of 
orbits. As expected, the Lyapunov exponents of orbits previously identified as quasiperiodic were found to vanish. 
The two orbits with an irregular broadband power spectrum, on the other hand, have finite and positive maximal 
Lyapunov exponents whose values were approximately determined as Lmax.i ^ 0.02 and L,nax,2 ~ 0.008. Those provide 
our most unequivocal and quantitative evidence for the chaoticity of the dyon-dyon interactions. The orbit with the 
smaller Lyapunov exponent shows in addition signs of intermittent behavior. 

Reassuringly, the results of all three employed analysis methods, i.e. Poincare sections, power spectra and maximal 
Lyapunov exponents, are fully consistent with each other. Taken together, they provide convincing evidence for and 
a quantitative description of both quasiperiodic and chaotic regions in the low-energy phase space of two BPS dyons. 
Moreover, the integrability of noninteracting dyon systems allows to trace the origin of the chaotic behavior to the 
interactions between the dyons. These interactions may therefore help to disorder monopole ensembles similar to 
those which are expected to populate the vacuum of the strong interactions. In any case, our results imply that no 
more than the three explicitly known integrals of the motion are conserved by the geodesic forces between the dyons. 
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APPENDIX A: ATIYAH-HITCHIN METRIC 



In this appendix we establish our notation and briefly summarize pertinent features of the Atiyah-Hitchin (AH) 
metric and its ChristofFel connection on the internal moduli space Mj"'' of the Bogomol'nyi-Prasad-Sommerficld (BPS) 
monopole pair. Full details can be found in the book [l^. 

BPS magnetic monopolcs are regular classical soliton solutions of the SU(2) Yang-Mills-Higgs equations in the 
limit of vanishing Higgs self-coupling [T^ IT3l |. The static two- monopole solutions define, as explained in Sec. |nl the 
cooUective-coordinate or moduli manifold Mj""* . Its metric was first written down explicitly by Atiyah and Hitchin js^ 
and determines the geodesic dynamics of two interacting BPS dyons at small velocities. The AH construction exploits 
the facts that the two- monopole moduli space is a hyper-Kahler (or Hamiltonian) manifold [Til Is^ , that it admits 
SO (3) as a group of isometrics and that the orbits under the action of SO (3) are with one exception three-dimensional. 

One can then show that the metric on A/j*'"' must be of the form 

ds^ = f^dg'^ + [a^lah + b^mamp + c^naUp) dx^'dx^, (Al) 

where a, b, c and / are fimctions of the "radial" variable g only and ladx"^, madx" and Uadx" are differential forms 
on the three-sphere which can be taken as 

(Ti I adx" = — sin t{.' d-d + sin d cos tp dif, 
a 2 '■= niadx"' = cos ip d'ff + sin 'd sin jJj dip, 

"■3 Uadx" = dip + cos^ dip, (A-2) 

in terms of three Euler angles 1}, ip and tp in the intervals 0<'i?<7r, 0<(^< 27r, < ip < 2tt. Accordingly, the 
two-monopole moduli space is parameterized by a coordinate g which describes the separation between the monopoles, 
two angular coordinates t9 and (p which determine the orientation of the axis that joins the monopoles, and the angle tp 
which fixes the position of the (generally axially asymmetric) two-monopole system with respect to rotations around 
this axis. 

By casting the line element ljAl|) into the Riemannian form gaf3 dx'^dx^ , one can straightforwardly verify that the 
fundamental metric tensor gap is symmetric, i.e., gap = gpa, and that the contravariant tensor g"'^ is its inverse, 
Qap g^^ = 5a'^ , as it should be. Explicitly, the diagonal elements of gap read 

900 "^f, 533 =c^, 
gii ~ c? sin^i/' + ^ cos^tp, 

922 = sin^?9 cos^-0 -|- sin^-d sin^7/i -I- cos^'d, (A3) 

while the nondiagonal ones are 

501 = 502 = .903 = 513 = 0, 

512 = {b'^ ~ a^) sin sin cos ■i/', 

523 = c^ cos I?. (A4) 
Similarly, for the elements of the contravariant tensor g"^ one has 

00 _ 1 33 1 , cosV , sin^V 



r-j2^ 5- = ^ + ^^ + ^)cot^^, 



5" 



cos^i/; sin^ip 
cos^ip sm^ip 



5'' = ( ^ + ^ ) csc^^, (A5) 



and 



g'' 



5" 



{c? — b^) CSC d sin tp cos ip 
^252 ' 

(b^ — a^) cot -d sin ip cos ip 



a^b'^ 

no ; cos^ip sin^iiN , . , 

3 =-(^^ + -^)csc?9coti?. (A6) 
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With these expressions at hand, one can obtain expUcit relations for the functions a, 6, c and / by means of the 
equation 



Act 



dx<^ 



i/a-L Act ~ i ua^ Xcx — U 

r(0) 



(A7) 



for the symmetric second-rank tensor R\cr which represents the fact that M2 is hypcr-Kahlcr. The V^x^ are the 
components of the Christoffel connection, i.e. the Christoffel symbols of the second kind, which are related to the 
metric tensor by 



\ dx^ dx"^ 



(A8) 



(Obviously, they are symmetric under exchange of the lower indices.) After inserting the expressions (|A3p - HA6p for 
the elements of the metric into Eq. (|A7|) . one obtains the components of R\„ as functions of a, 5, c and /: 



a' 


b' 




a" 


b" 


c" 


a 




4)f 


a 


b 


c 



Roa 

(the prime denotes differentiation with respect to the coordinate g) and 

Rii = Ui sin^V + n2 cos^Vj 
R12 = (112 — Hi) siniJsin'i/'cos'!/', 
R22 = (ill cos^i' + II2 sin^ip) sin^z? + U3 cos^t?, 
^33 = Ha, 



(A9) 



R 



23 



n.s COS-&, 



(AlO) 



where 



Hi 



n. 





262c2 












2a? 




c4 







2a262 



a 



b' 



c — c 



a 
b 

c 

P' 



(All) 



All other components of R (with the exception of those which differ from the above by exchanging the indices) vanish. 
The equations (|A7|) therefore reduce to 



Hi =0, Ha = 0, = 0, 



together with 



al_b^ 
a b 



'b~ 



1 

62 



20252^2 



(A12) 



(A13) 



The last equation, however, is just a first integral of Eqs. ljA12ll and may be regarded as a constraint imposed on the 
initial values of a, 6, c and its derivatives. The second-order equations, which of course conserve this constraint, can 
be obtained from each other by cyclic permutation of (a, b, c). All four equations ljA12p and IjAlSp together constitute 
the vacuum Einstein equations of the AH metric. 

A particular set of first integrals of the second-order equations l)A12|l are the first-order equation 



2bcda 2 



(AM) 



and the two others which are obtained from it by cwlical permutation of (a, 6, c). These differential equations for the 
functions a, b, c and / were first derived in Ref. (57|- They have been linearized and solved in terms of Legcndre's 
complete elliptic integrals of the first and second kind in Ref. [s^ for / ~ abc and in Ref. for / = —b/g. Below 



13 



we will adopt the second choice, / = —b/g, since it leads to expressions which arc more convenient for our purposes. 
The explicit solution to Eq. (|A14() and its permutations is then 

= U\{K -E){E~- q^K)/E, 
b^ = 4:KE{K - E)/{E - q^K), 

= AKE{E - fK)/{K - E), (A15) 

where 

fw/2 i-7r/2 

K{q)= dT{l~ q^sm^T)-^''^, E{q) = dr (1 - q^ sin^r)^^^ (A16) 

Jo Jo 

are the complete elliptic integrals of the first and second kind, q is related to g hy g = 2K(q) for tt < g < oo and 
q = — is the conjugate modulus. At the value g = it (which implies q = 0) the metric has a coordinate 
singularity, the so-called "bolt", since a(7r) = implies that the line element IjAip becomes independent of cti. 



APPENDIX B: EQUATIONS AND INTEGRALS OF MOTION 
1. Lagrange equations of motion 

The geodesic low-energy dynamics of the two-dyon system is governed by the Lagrangian 

L = jgo.,i"i^ - I { f(^' + aW, + bW^ + cW^) , (Bl) 

where the reduced mass m of the dyon-pair has been set to to = 1 (instead of to = 27r in Ref. The first equation 

in (|Bip describes generic geodesic motion while the second one is specialized to the Atiyah-Hitchin metric on the 
internal collective coordinate manifold Mj"'' in terms of the functions a, b, c and /. The components lo^, ujy and 
of the instantaneous (or "body- fixed" ) angular velocity Co along the axes x, y and z may be expressed in terms of the 
rates of change of the Euler angles as 

LOx = 'Ji/dt^ LUy = a2/dt, LOz — a^/dt (B2) 
(cf. Eq. (|A2|) '). The canonically conjugate momenta associated with the four relative coordinates g, i?, (p and ip are 

Pg = f^Q, Pi' = 

Pf) = b^ojy costA — cL^^jJx sinV', 

Pip ~ a^ujx sin ?9 cos -t- b'^bjy sin sin + (?uJz cos d. (B3) 

The Euler-Lagrange equations for the time evolution of the collective coordinates are obtained by varying the action 
based on the Lagrangian (|ljl|) . Variation with respect to the radial coordinate g{t) leads to 

fQ + ffo" - a a'iol - b b'to^ - c c'lu^ = 0, (B4) 

while the equation of motion for the coordinate ijj{t) becomes 

2cc'iLJzg + c^i^z — {a^"^ — b~'^^ a^b^LOx^y (B5) 

The remaining two equations of motion for Lp{t) and i9(t) are 

2b b' Ldyi) cos ^jj — 2a a'ujxgsin^ + c^uJz^sind 
= (iJjx sin ip + ijJx^z cos -0) a? — {6jy cos — uiytUz sin -0) b^ (B6) 

and 

2 a a'uJx g sin cos ip + 2b b'ujy g sin 1) sin -I- 2 c c'ujz g cos 
= [ (0 sin 1? sin -0 — ?9 cos i? cos -0) uJx — dJx sin i? cos -0] 
— [(4! sin I? cos -0 -f- 1? cos i? sin ^/j) ojy + ujy sin i3 sin -0] b^ 

+ (ujz^ sin ^ -ojz cos 1?) . (B7) 
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A pair of equations which are algebraic consequences of Eqs. (|B6|) and (|B7|I . 

2bh'ujyQ + h'^L0y = (c^^ - a"^) a^c^w^jW^, (B8) 
2aa'ujxQ + a^tja- = {b^^ — c^^) b^c^ujyUJz, (B9) 

will be helpful in arriving at a more efficient formulation. Indeed, with q = Pg/ f'^ one also has 



and, introducing the (^?-dependently) rescaled angular velocities (in the body-fixed frame) 

Ml = a^uj^, M2 6^tJy, M3 = c^w^, (Bll) 
one can rewrite the four equations of (relative) motion concisely as 

Ml = (^^ - M2M3, (B12) 

M2- (^^-^^MaMi, (B13) 

M3- (^^-p^MiM2, (B14) 

P, = ^ + M?^+M^^+Mi^. (B15) 
In this form, the equations of motion were first obtained in Refs. pflllH^. 



2. Integrals of the motion 



Three independent [821 constants of the geodesic Atiyah-Hitchin motion are known explicitly. The first one can be 
immediately identified by noting that (/? is a cyclic coordinate, i.e. that it does not appear explicitly in the Lagrangian 
(|Blll . The corresponding generalized momentum p^p is therefore an integral of the motion, i.e. p^p ~ 0. 

The Hamiltonian H is obtained by the standard Legendre transformation 



of the Lagrangian (jBl|l . Using Eqs. 
momenta, the Hamiltonian becomes 



H = pgQ + p^'d + p^ip + p^,^j - L (B16) 
to express velocities in terms of coordinates and their conjugate generalized 



H 



2 



„2 



sin^-i/j 



cos^i/; 



COS^l/' 



^^^^ ) (p^ CSC ■& - p^ cot z?)^ 



Pi) 



{pip CSC d ~ p^ cot -d) sin (2?/') 



(B17) 



which assumes the same values as the Lagrangian since both are of purely kinetic origin. Neither in L nor, as 
a consequence, in H does the time coordinate t appear explicitly. The corresponding Hamilton equation, H = 
{H,H} = 0, therefore trivially implies that the total relative energy H of the two-dyon system is conserved and 
furnishes a second integral of the motion. 

The third constant of the motion is the square of the (both frame- and body-fixed) total angular momentum, M^ , 
which can be expressed as the sum of the squares of the rescaled body-fixed angular velocities, i.e. M^ = Mf -|-M2-|-M|. 
One may check that M^ is conserved by adding up the equations obtained from multiplying Eq. HB12|I by 2Mi, 
Eq. ||BT3)| by 2M2 and Eq. l|BT4jl by 2M3. The result is 



2M1M1 + 2M2M2 + 2M3M3 = 



and therefore 



dt 



M^ = 0. 



(BIS) 
(B19) 
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By using the explicit expressions 



PipCosV'-pvCOSi?cos'0-p^sinz?sinV' 

JYLi = : — , [tiZU) 

sm w 



Pip sin ip ~ P4> cos 1? sin + sin cos -i/; 

sin-d ' 



(B21) 



and 



M3 = p^ (B22) 

for the rcscalcd angular velocities, one can re-express the total angular momentum squared in terms of the Euler 
angles and their canonically conjugate momenta as 

3 

= ^ M| = - 2p^p^ cot 7? CSC 7? + {pI + pI) csc^d. (B23) 

For the investigation of chaos in the two-dyon system it is important to note that a fourth integral of the motion 
appears in the g —^ 00 limit, i.e. if the dyons are far separated. Indeed, in this case the Atiyah-Hitchin metric 
simplifies to the Euclidean Taub-NUT metric and the electric charge density, associated with the canonical 
momentum p^ = c^ujz, becomes a fourth constant of the motion, i.e. p^ = 0. 

The simple physical explanation of this result is that infinitely separated dyons cannot exchange electric charge, 
so that in addition to their overall charge also their individual charges become time- independent. In this case there 
exist at least four independently conserved quantities |83l | in the eight-dimensional phase space, enough to render 
the geodesic dynamics in Euclidean Taub-NUT space (Liouville-) integrable As a consequence, the asymptotic 
motion of two BPS dyons cannot be chaotic. Moreover, the potential transition from regular motion at g — > 00 to 
chaotic motion for finite, decreasing g may be "delayed" according to KAM theory [sj by a stepwise dissolution of 
the invariant tori. 
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[61] The chaoticity of gauge field trajectories was characterized by positive maximal Lyapunov exponents 0, the whole Lya- 
punov spectrum 0, and the Kolmogorov-Sinai entropy |3- In combination, these properties indicate the global hyper- 
bolicity of classical non-Abelian lattice gauge theory (in 3+1 dimensions). Other important findings were a continuous 
cascading of the dynamical degrees of freedom (and their energy) towards the ultraviolet during time evolution 0, which 
might impede the continuum limit of non-Abelian lattice gauge theories j^, and the distribution of nearest-neighbor level 



17 



spacings of Yang-Mills and QCD lattice Dirac spectra according to Gaussian matrix ensembles , indicating from another 
perspective that the underlying classical theory is chaotic 

[62] At high temperatures, long-wavelength fields behave increasingly classical. Chaos investigations can therefore provide useful 
input to the study of otherwise hardly accessible nonequilibrium processes in hot gauge theories. Timely applications include 
calculations of the fast thermalization rates observed in heavy-ion collisions jlj] and (especially topological) structure 
formation, e.g., in baryon number violating processes during semi-classical evolution phases in the standard model. 

[63] Prominent examples include the interactions of Skyrmions, vortices, instantons and D-branes. 

[64] As a case in point, in the 2-1-1 dimensional Yang-Mills-Higgs model 't Hooft-Polyakov monopoles (a generalization of 
BPS monopoles which play the role of instantons in this case) generate "weak confinement" by forming a monopole 
antimonopole plasma, as shown long ago by Polyakov [l9l[ . and the "strong confinement" of 2+1 dimensional Yang-Mills 
theory is expected to be due to a similar mechanism. 

[65] Such solutions may be similar or related to the BPS dyon constituents of the recently discovered SU(N) Yang-Mills caloron 
solutions with nontrivial holonomy [25II . (In the case of SU(2), they contain a BPS monopole-antimonopole pair.) 

[66] At fixed time, the polar angles i9 and specify the direction of the axis connecting the two monopoles, and t/j corresponds to 
rotations around this azxis. The angle if) becomes relevant since the BPS monopole pair is in general not axially symmetric. 

[67] despite the long-range forces generated by the massless Higgs fields 

[68] Recall that this motion is not ergodic since the remaining three constants of motion prevent it from filling the seven- 
dimensional constant-energy surface densely and uniformely. 

[69] The existence of bounded quasiperiodic geodesies in the almost asymptotic AH metric, i.e. the existence of quasiperiodic 
solutions to the dynamics ^ in its sufficiently weakly nonintegrable realms, has been established in Ref. l3al and shown 
to form a set with positive Lebesgue measure. 

[70] More specifically, a simple periodic motion leads to a single fixed point in the Poincare section while a periodic orbit 
with two commensurable frequencies (whose quotient is a rational number) gives rise to a finite number of points being 
indefinitely repeated in the same order. A two-frequency quasiperiodic orbit draws a closed curve that never exactly retraces 
its steps, and chaotic motion appear as a scatter of points which eventually fill a two-dimensional surface. 

[71] In the large-g region, the dyon-pair orbits can be obtained analytically by perturbation of the asymptotic Taub-NUT 
metric and turn out to be conic sections [l5l| . It would therefore be possible to obtain analytical expressions for the 
Poincare sections in this region, too. Since we expect the chaotic regime to develop at smaller g values, i.e. outside the 
range of validity of the (first-order) analytical results, however, we will resort to a numerical treatment of all orbits. 

[72] Note that the times between consecutive piercings of the Mi = plane by an orbit typically vary. This is in contrast to 
stroboscopic studies where one determines the system's output at equal time intervals. 

[73] Numerically, this can be checked by calculating the time dependence of along the trajectory. For the chaotic orbit it is 
indeed much stronger then for all previous ones. 

[74] Note that we plot the logarithm of the power spectra in order to emphasize the structure of the background. 

[75] Rational independence of v\,2 implies that the only solution of miV\ -\- m2V2 = is mi = 1x12 = 0. 

[76] The choice of the base frequencies is to a certain extent a matter of convention |4^ . 

[77] We recall the well-known effect that in a continuous power spectrum based on a finite time interval the contributions from 
the lowest frequencies are artificially enhanced. This is due to the fact that aperiodic points in a finite data set appear as 
points with very long periods, comparable to the length of time over which the series is recorded. 

[78] In fact, it appears that the quasiregular orbits unfold the skeleton of the aperiodic signals. A further analysis of such orbits 
could be persued by the method of Ref. li^l which allows to separate the power spectra of conservative systems into sharp 
and broadband features by utilizing the distribution of local Lyapunov exponents lisll . 

[79] It would be interesting to study the onset of chaos further, e.g., by numerically identifying orbits on commensurate- 
frequency KAM tori which nonintegrable perturbations destroy first. 

[80] This holds even more strongly for dissipative systems with strange attractors because power spectra discard the phase 
information of the Fourier spectrum (cf. Eq. 11131 '). 

[81] For a more detailed understanding of chaotic systems, the remaining (nonmaximal) Lyapunov exponents are also of interest. 
Their numerical evaluation is significantly more involved, however lEH. l52l Is^l . 

[82] Two integrals of the motion are independent if their mutual Poisson bracket vanishes. 

[83] In fact, there are even more constants of the motion j^. Those are associated with the S0(4) (SO(3, 1)) symmetry of 
bounded (unbounde d) g eodesic motion in Euclidean Taub-NUT , which can be traced to the existence of a Killing- Yano 
tensor in this space 1601 . 
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FIG. 2: Same as Fig.Elwith = 2152.95. 




FIG. 3: Same as Fig. 



0with = 2237.7. (The outermost section contains about 13500 points.) 
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FIG. 4: Panel a (b, c) shows the time dependence of pg and panel e (f, g) the logarithm of the corresponding power spectrum 
for the dyon-pair orbit whose Poincare section is plotted in Fig.Ql|21El- The panels d and h show Pgit) and the corresponding 
logarithmic power spectrum for the additional chaotic orbit described in the text. 




FIG. 5: The full (dotted, dash-dotted, dashed) curve depicts x(^) for the dyon-pair orbit whose power spectrum is plotted in 
Fig. He (f, g, h). 
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FIG. 6: The function for three quasiperiodic dyon-pair orbits. The full (dotted, dash-dotted) curve corresponds to the 
orbit with the innermost Poincare section in Fig. . 



